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KMAX 
N 1 
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N4 
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SX1 

SX2 

SX3 

SYO 

SY1 

SZO 


NOMENCLATURE 


length  of  afterbody 

distance  from  leading  edge  of  airfoil  to  initial  value  line 

distance  from  tail  of  body  to  outflow  boundary 

distance  from  trailing  edge  of  airfoil  to  tail  of  afterbody 

number  of  grid  points  on  stagnation  line 

number  of  stacked  surfaces,  radial  direction 

number  of  grid  intervals  on  airfoil  between  leading  edge  and 
point  K 

number  of  grid  intervals  on  airfoil  between  point  K  and  trailing 
edge 

number  of  grid  intervals  on  afterbody  between  airfoil  trailing 
edge  and  end  of  body 

number  of  grid  intervals  between  tail  of  body  and  outflow 
boundary 

grid  clustering  parameter  at  airfoil  leading  edge 
grid  clustering  parameter  at  point  K  on  Segment  1 
grid  clustering  parameter  at  airfoil  trailing  edge 
grid  clustering  parameter  at  tail  of  body 

grid  clustering  parameter  on  stagnation  line  at  airfoil  leading 
edge 

grid  clustering  parameter  on  stagnation  line  at  initial  value 
surface 

grid  clustering  parameter  on  afterbody  surface  in  radial 
direct  ion 


All  other  symbols  are  defined  in  the  text 
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1 .  Introduction 

Algebraic  grid  generation  nethods  for  three-dimensional  (3-D)  flow 
problems  have  the  advantage  over  their  differential  equation  counterparts  in 
speed  and  ability  to  handle  high  aspect  ratio  cells  without  difficulty. 

Where  the  algebraic  methods  are  sometimes  at  a  disadvantage  is  in  treating 
a  wide  variety  of  boundary  shapes  with  a  single  code. 

Since  1979  grids  about  wing-body  configurations  have  been  successfully 
generated  by  several  algebraic  approaches.  Eriksson  [1]  has  generated  a 
single-block  nonorthogonal  3-D  grid  using  transfinite  interpolation  where 
geometric  data  is  specified  only  on  the  boundaries.  Since  no  internal 
surfaces  are  specified,  grid  quality  is  controlled,  especially  near  a 
surface,  by  incorporating  out-of-surface  parametric  derivatives.  Smith  [2] 
uses  the  patched  grid  approach  where  the  domain  is  divided  into  regions 
with  boundaries  of  a  simpler  character  than  the  overall  region.  On  the 
interior  of  each  six-sided  sub-region  transfinite  interpolation  is  used  to 
generate  the  grid.  His  treatment  extends  only  to  the  wing  tips  which  limits 
its  usefulness.  A  third  and  quite  different  approach  has  been  taken  by 
Caughey  &  Jameson  [3].  Their  technique  generates  a  boundary-conforming 
coordinate  system  by  a  sequence  of  conformal  and  shearing  transformations 
to  yield  a  nearly  orthogonal  computational  domain.  The  grid  is  then 
generated  by  simple  linear  interpolation.  Shmilovich  &  Caughey  [4]  have 
recently  extended  this  technique  to  include  a  tail  surface.  The  Caughey- 
Jameson  procedure  was  developed  for  use  with  the  3-D  transonic  FhO  codes. 

One  of  the  maior  difficulties  in  algebraic  grid  generation  is  preventing 
corner  singularities  on  the  boundaries  from  propagating  into  the  grid.  Any 
i nterpolation  method  will  propagate  such  singularities  into  the  interior. 
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Differential  equation  grid  generation  schemes  suffer  no  such  problem  because 
of  the  diffusive  action  of  the  elliptic  operator.  Corner  singularities 
are  always  present  in  3-D  configurations  at  wing/fin  trailing  edges  and 
tail  points  on  pointed  bodies.  A  method  for  removing  these  singularities 
in  algebraic  grid  generation  has  been  developed  by  Vinokur  &  Lombard  [5]  for 
2-D  geometries.  Their  method  consists  of  patching  a  conformal  hinge  point 
transformation  in  a  small  region  near  the  corner  to  a  grid  in  the  remainder 
of  the  domain  generated  by  transfinite  interpolation.  They  successfully 
applied  this  method  in  generating  a  patched  grid  in  a  domain  consisting  of 
a  backward  facing  step  at  the  end  of  a  nozzle  exhausting  into  a  cylindrical 
diffuser . 

The  problem  treated  in  this  report  is  the  generation  of  a  surface  fitted 
grid  in  the  stern  region  of  an  undersea  vehicle,  specifically  an  axisymmetric 
pointed  afterbody  with  four  identical,  symmetric,  constant  chord  fins.  In 
many  respects  this  problem  is  similar  to  the  wing-fuselage  problem.  The 
desired  grid  is  to  be  used  for  either  inviscid  or  viscid  incompressible  flow 
calculations  and  hence  must  have  proper  clustering  ability  to  resolve 
regions  of  high  flow  gradients.  An  algebraic  approach  is  used  which  is  an 
outgrowth  of  earlier  3-D  grid  generation  work  on  a  fin-cylinder  configuration 
[6]  . 

The  algebraic  method  adopted  here  was  originally  inspired  by  the  work  of 
Caughey  .5  Jameson  in  unwrapping  a  geometrv  as  much  as  possible  to  produce  a 
naral 1 elpiped  with  nearly  straight  boundaries.  This  procedure  is  ot  the 
stacking  type  where  a  3-D  grid  is  produced  by  a  sequence  of  2-D  grid 
generation  operations.  in  the  present  method  stacked  tubular  surfaces  ot 
circular  cross-section  are  first  determined  and  then  a  C-type  grid  generated 
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on  each  surface.  In  the  process  of  generating  these  surfaces  as  well  as 
in  unwrapping  the  airfoil,  corner  singularities  are  removed  by  application 
of  a  hinge  point  transformation  to  the  entire  boundary.  The  present 
approach  thus  differs  from  that  of  Vinokur  &  Lombard  in  being  global  rather 
than  local  in  the  use  of  the  hinge  point  transformation.  The  result  is  a 
smooth  boundary  with  a  slowly  varying  tangent.  A  grid  which  is  orthogonal 
at  all  boundaries  is  then  generated  on  the  interior  by  transfinite 
interpolation.  By  using  a  sequence  of  one-dimensional  stretching  functions 
in  physical  space,  precise  control  is  maintained  over  the  clustering  at 
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2 .  Analysis 

2 . 1  Geometry  of  Computational  Domain 

We  start  by  defining  the  geometry  for  which  a  surface-fitted  grid  is  to 
be  generated. 

(1)  The  afterbody  is  of  circular  cross-section  and  has  a  smooth 
but  otherwise  arbitrary  meridian  profile  that  closes  at  the 
tail  point. 

(2)  Four  identical  fins  of  constant  unit  chord  and  infinite 
span,  consisting  of  symmetric  airfoil  sections,  are 
mounted  90  degrees  apart  with  their  chord  planes  passing 
through  the  afterbody  centerline.  The  trailing  edges  of 
the  fins  are  located  upstream  of  the  tail  point  a  distance 
dTL* 

(3)  The  computational  domain  consists  of  the  region  interior  to 
an  outer  cylinder  of  radius  r^xp  and  exterior  to  the 
afterbody,  bounded  upstream  and  downstream  by  planes  normal 
to  the  afterbody  centerline  (the  initial  value  and  outflow 
planes) . 

A  schematic  of  the  geometry  and  computational  domain  is  shown  in  Fig.  1 
and  a  head-on  view  showing  the  coordinate  system  in  the  crossflow  plane 
appears  in  Fig.  2.  Since  the  fins  are  identical  and  equally  spaced  there 
are  four  planes  of  symmetry  at  9  =  0,  n/4,  tt/2  and  3tt/4.  Thus,  only  the 
section  -  tt/4  <  9  <  0  is  considered  in  generating  a  grid  and  in  the  flowfield 


calculation. 
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2.2  Grid  Stacking  Procedure 

The  simplest  grid  stacking  scheme,  such  as  that  in  Ref.  3,  makes  use  of 
a  shearing  transformation  to  distribute  a  sequence  of  two-dimensionally 
produced  grids  in  the  third  spatial  direction.  Unfortunately  the  shearing 
transformation  causes  surface  corner  discontinuities  to  propagate  into 
the  grid.  In  the  present  case  of  a  pointed  tail  body,  Fig.  3  illustrates 
the  situation  that  would  exist  in  the  meridian  plane  if  a  shearing 
transformation  were  used.  Along  the  vertical  line  through  the  tail  point, 
x  =  x-p,  lines  of  constant  n  have  discontinuous  slopes. 

What  is  needed  is  a  transformation  to  produce  n  =  constant  lines  that 
does  not  propagate  corner  discontinuities.  The  hinge  point  (power  law) 
conformal  transformation  has  this  desired  property.  At  the  tail  point 
(corner  discontinuity)  we  write 


where  the  real  axis  is  aligned  with  the  axis  of  symmetry,  and 

z  =  x  +  iy  , 

w  =  u  +  iv  , 

TT 


(1) 

(2) 

(3) 

(4) 


A 

The  taiL  angle  0^,  is  defined  in  the  meridian  plane  schematic  given  in  Fig.  4 
Equation  (1)  maps  the  sector  0  <,  0  <  n  -  9  above  the  real  axis  onto  the 
upper  half  pLane.  If  Eq .  (1)  is  applied  to  the  entire  bounding  curve  in 
the  meridian  plane,  A-B-C-D-E-A,  the  corner  at  the  tail  point  B  is 
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eliminated  in  the  w  (hinge)  plane.  Then  interpolating  a  grid  in  the  hinge 
plane,  upon  transformation,  will  produce  a  smooth  grid  (except  at  B)  in  the 
physical  plane  which  can  be  used  for  stacking. 

Before  proceeding  further,  we  present  the  calculational  steps  in  going 
from  the  z  plane  to  the  w  plane  and  vice-versa.  We  start  by  writing  z  and 
w  in  polar  form, 


z 


-  16  14 

re  ,  w  =  pe 


(5) 


Then  substituting  Eqs.  (5)  into  Eq.  (1)  and  equating  real  and  imaginary  part: 
gives 


P 


«n 

r 


$  = 


a 

n9 


> 

J 


(6) 


Given  (x,y),  to  compute  (u,v)  we  first  compute  the  magnitude  r  and  angle  0 
from 


A 

0  = 


i 


TT 

2 


,  x  =  0 

,  x  >  0 


(7) 


Next,  p  and  $  are  computed  from  Eq.  (8)  and  finally  u  and  v  (transforming 
from  polar  to  cartesian  form)  from 
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►  k  ' 


m 


1  V. 


u  =  p  cos 


v  =  p  sin  $ 


(8) 


Conversely,  if  (u,v)  are  given,  the  magnitude  p  and  angle  <fi  are  computed 


from 


r  2  211/2 

p  =  [u  +  v  J 


tan  *(— 1  +  tt  ,  u  <  0 


IT 

2 


tan  1  (— ) 
v  u' 


,  u  =  0 

,  u  >  0 


Then  r  and  0  are  computed  from  the  inverse  of  Eq.  (6),  viz. 


r  =  p 


l/n 


6  =  —  4> 
n 


Finally,  x  and  y  are  related  to  r  and  0  by 


x  «*  r  cos  0 


v  =  r  sin  0 


The  image  of  boundary  A-B-C-D-E-A  in  the  hinge  plane  is  shown 


(9) 


(10) 


(ID 


schematically  in  Fig.  5  and  has  the  appearance  of  a  water  spout.  Segment  BC 
remains  straight  because  the  real  axis  in  the  z  plane  is  coincident  with 
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BC.  For  convenience,  the  transformed  boundary  segments  will  be  denoted  hv 
noint-of-compass  notation.  Thus  A-B-C  is  v,,(u)  and  E-D  is  v^(u). 

To  accommodate  the  inflow  and  outflow  boundaries  (lines  A-E  and  C-D  in 
Fig.  5)  as  well  as  the  leading  and  trailing  edges  of  the  fin,  vertical  lines 
must  remain  vertical  in  the  physical  plane.  This  requirement  acts  as  a 
constraint  on  the  transformation  from  the  hinge  to  the  physical  plane. 

The  simplest  scheme  for  producing  a  grid  in  the  hinge  plane  is  a  shearing 
transformation  on  the  image  of  x  =  constant  lines.  Thus  the  normalized 
variable  r  is  defined  as 


x=const . 


(12) 


At  this  point  the  distribution  of  r.  is  assumed  known.  Thus  in  the  interior 

J 

v  is  given  by 


v 


+ 


9 


(13) 


where  the  index  i 
Thus  at  point 
are  determined  by 
for  y  is  obtained 


is  constant  on  x  =  constant  lines. 

(i,j)  the  values  of  x.  and  v.  .  are  known.  Then  (y,u).  . 

l  i,j  J  i,j 

iteratively  solving  Eq.  (1)  as  follows:  A  first  guess 
by  linear  interpolation  from 


VyN 


(14) 


The  cycle  begins  by  computing  u  from 
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A  A 

where  r  and  9  are  given  by  Eqs.  (7).  Next,  new  values  of  p  and  <j>, 
corresponding  to  the  current  estimate  of  u  and  the  known  v,  are  computed 
from  Eqs.  (9).  The  p  and  4>  are  then  used  to  update  r  and  9  using  Eqs.  (10), 
and  finally,  xc  and  y  are  given  by  Eqs.  (11).  Convergence  is  attained  when 

I x  -  xc|  <  10  6  .  (16) 

Thus  for  each  value  of  r  =  constant,  a  smooth  curve  y^(x)  is  determined 
in  the  meridian  plane  using  the  above  procedure.  By  revolving  y^,(x)  about 
the  x-axis  a  tubular  coordinate  surface  is  obtained  which  is  smooth  and  non- 
developable  (except  when  it  is  a  cylinder).  On  each  of  these  surfaces  a 
surface  fitted  grid  is  determined  as  though  the  surface  were  developable, 
then  projected  back  onto  the  surface.  This  means  that  given  (x,6),  r  is 
determined  by  interpolation  of  the  tubular  meridian-plane  curve  rg(x)  =  yg(x). 
For  the  projection  method  to  work  properly,  the  foil  subtended  angle  9  (x) 

f 

must  be  computed  to  account  for  the  variable  r<,(x)  from 

y 

0F  -  sin-1(/)  (17) 

S 

where  yp(x)  is  the  airfoil  semi-thickness  distribution.  Lagrange  cubic 
interpolation  is  used  to  determine  r,  given  x,  from  the  previously  determined 
values  of  r^. 

Clustering  of  r  =  constant  lines  near  the  body  surface  is  needed  to 
resolve  the  viscous  layer  whereas  further  away,  where  flow  gradients 
diminish,  these  lines  can  be  further  apart.  A  one-sided  stretching  function 
is  therefore  appropriate  to  determine  the  grid  line  spacing  in  the  meridian 


plane 
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Vinokur  [7]  has  determined  approximate  criteria  for  the  development  of 
one-  and  two-sided  stretching  functions  of  one  variable  which  give  a  uniform 
truncation  error  independent  of  the  governing  differential  equation  or 
difference  algorithm.  He  investigates  several  analytic  functions  but  finds 
that  only  tan  z,  where  z  is  real  or  pure  imaginary,  satisfies  all  of  his 
criteria.  These  stretching  functions  were  used  in  the  predecessor  grid 
generation  scheme  [6]  and  are  also  used  here. 

Since  r  is  already  a  normalized  variable,  its  distribution  is  given  by 


tanh[y  -  l)] 

tanh 


(18) 


where  £  is  the  normalized  generating  variable  given  by 


and  A<j>  is  the  solution  of 


S0  = 


sinh  A(j> 
A<{> 


and 


*. '  7 <0> 

dr 


(19) 


N-  =  number  of  intervals  in  r  =  KMAX  -  L 


2 .3  Grid  Generation  on  a  Tubular  Surface 
Overall  Approach 

Grid  generation  in  the  x  -  9  plane  is  accomplished  in  three  stages.  The 
first  stage  involves  a  sequence  of  conformal  transformations  to  unwrap  the 
airfoil,  symmetry  lines  and  initial  and  outflow  lines  into  a  quadrilateral 
with  a  slowly  varying  height.  The  unwrapping  transformations  are  the  basis 
for  producing  a  C-grid  about  the  airfoil.  The  second  stage  involves 
translation  and  rotation  of  coordinates  about  the  image  of  the  airfoil 
trailing  edge,  followed  by  a  hinge  point  transformation  to  eliminate  the 
corner  at  the  trailing  edge.  The  third  stage  makes  use  of  transfinite 
interpolation  to  determine  the  grid  in  the  hinge  plane  that  is 
orthogonal  at  all  boundaries.  Since  the  boundaries  in  the  hinge  plane  are 
smooth  and  have  a  slowly  varying  tangent,  transfinite  interpolation  will 
produce  a  smooth  grid  in  which  non-orthogonality  in  the  interior  is  held  to 
a  minimum.  The  grid  in  the  physical  plane  is  obtained  by  taking  the  inverse 
of  the  sequence  of  transformations.  Since  the  intermediate  transformations 
are  conformal,  the  orthogonality  at  boundaries  and  grid  smoothness  will  be 
preserved  in  the  physical  plane.  Spacing  of  grid  lines  is  determined  on 
appropriate  boundaries  in  the  physical  plane  by  use  of  stretching  functions. 

Sequence  of  Transformations 

At  this  stage  the  coordinates  of  the  fin  (x  ,9  )  on  a  specified 
r  =  constant  tubular  surface  are  assumed  known.  The  first  step,  in 
preparation  for  the  unwrapping  transformation,  is  to  scale  the  (x,0) 
coordinates  according  to 
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x  =  4(x  -  d  )  +  in  2 


0  =  49 


(20) 


where  dg  is  the  location  of  the  singular  point  in  the  unwrapping  trans¬ 
formation  and  is  just  inside  the  leading  edge  of  the  airfoil.  The  stretching 
factor  4  is  required  by  the  unwrapping  transformation  so  that  the  upper  limit 
on  0  will  be  ±  it  (the  upper  and  lower  symmetry  planes  are  at  0  =  ±  ir/4). 

In  Refs.  3  and  6,  x  is  translated  but  not  magnified  whereas  0  is  magnified 
as  above.  Unequal  scaling  is  of  course  not  conformal  so  that  orthogonality 
of  the  grid  cannot  be  maintained  at  the  boundaries.  The  resulting  grid  in 
the  x  -  0  plane  will  be  highly  flattened  and  thus  highly  nonorthogonal . 

On  an  r  =  constant  surface  the  boundaries  and  coordinate  system  in  the 
x  -  0  plane  are  sketched  in  Fig.  6.  Because  of  symmetry,  only  the  region 
-  it  <  9  <  0  needs  to  be  considered.  The  airfoil  can  be  unwrapped  by  applying 
the  conformal  transformation, 

x  +  i9  =  £n[l  -  cosh(£  +  ip)]  .  (21) 

Equation  (21)  maps  the  region  below  the  x-axis  to  positive  K  in  the  band 
0  <  n  <  tt .  We  note  Chat  the  sign  of  9  in  Eq.  (21)  follows  Caughey  &  Jameson 
[3]  so  that  a  right-handed  coordinate  system  results.  In  Ref.  b,  the  sign  of 
9  was  negative  which  produced  a  left-handed  system.  Figure  7  presents  a 
schematic  of  the  boundaries  in  the  £  -  n  plane.  The  initial  value  line  (IVL) 


A-B-C  is  seen  to  map  into  a  near  semi-circle. 
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Following  Ref.  6,  the  corners  at  points  A  and  C  can  be  eliminated  by 
applying  the  conformal  transformation, 

-  -  C0 

5  +  in  =  £  +  in  +  -g-  ,  (22) 

where  £^  is  defined  in  Fig.  7.  Equation  (22)  has  the  effect  of  nearly 
straightening  out  the  IVL.  The  geometry  of  the  boundaries  in  the  £  -  q 
plane  is  shown  in  Fig.  8. 

In  Ref.  6  a  shearing  transformation  is  used  to  produce  a  grid  in  the 

£  -  n  plane.  For  airfoils  with  non-zero  trailing  edge  angles  this  procedure 

produces  discontinuous  metric  coefficients  across  the  line  £  =  £„.  To 

F 

eliminate  the  effect  of  the  corner  at  the  trailing  edge  (point  F)  a  procedure 
similar  to  the  generation  of  the  smooth  curves  in  the  meridian  plane  is  used. 
The  £  -  n  coordinates  are  first  translated  and  rotated  about  point  F  according 
to 


x  =  (£p  -  £)cos  Ap  +  (np  -  n)sin  Ap 
y  =  -  (IF  -  I)sin  Ap  +  (nF  -  h)cos  Ap 


(23) 


where  the  positive  x  axis  points  toward  point  E,  the  airfoil  leading  edge, 
and  Ap  is  the  trailing  edge  angle  in  the  £  -  n  plane.  The  translated  and 
rotated  £  -  n  boundaries  are  sketched  in  Fig.  9. 
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The  final  step  in  producing  a  smooth  boundary  is  to  apply  a  hinge  point 
transformation  to  (x,y)  to  eliminate  the  corner  at  point  F.  Equation  (1) 

A  A 

applies  provided  z  =  x  +  iy  and 


where 


n 


7T 

ir  -  AX  ’ 


(24) 


= 


(25) 


The  resulting  boundary  in  the  hinge  plane  is  sketched  in  Fig.  10.  Once  a 
grid  in  the  hinge  plane  is  produced  by  transfinite  interpolation,  the 
transformation  sequence  is  reversed  to  obtain  the  grid  in  the  x  -  9  plane. 

The  FORTRAN  coding  is  written  in  terms  of  real  variables  which  involves 
determining  the  real  and  imaginary  parts  of  the  conformal  transformations 
as  well  as  their  inverses.  These  relations  are  given  in  the  appendix. 


Corner  Point  Projection  on  Airfoil 

Part  of  the  present  grid  generation  strategy  is  to  force  one  of  the 
coordinate  lines  normal  to  the  airfoil  to  pass  through  point  C,  the  corner 
point.  This  point  on  the  airfoil  is  denoted  by  letter  K  —  see  Fig.  10. 
Such  a  line  provides  a  natural  division  between  those  lines  intersecting 
the  airfoil  from  the  IVL,  B-C,  and  the  lower  symmetry  line,  C-.J . 

An  effective  method  of  locating  point  K  that  prevents  reflexes  on  the 
connecting  segment  C-K  is  to  construct  a  circular  arc  between  C  and  K  which 
is  normal  to  both  boundaries.  Under  this  assumption  the  relation  between 
the  coordinates  at  C  and  K  is  found  to  be 


tan[i  f>c  +  4>k)  ] 


(2b) 


where 
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(27) 


Supplementing  Eq .  (26)  by  the  equation  for  the  airfoil  image,  v^(u),  gives 

two  equations  for  the  two  unknowns  u  and  v  . 

K  K 

Equation  (26),  together  with  the  airfoil  image  equation,  can  be  solved 
iteratively  by  the  following  formula,  derivable  from  Newton's  method: 


where 


The  right-hand-side  of  Eq .  (29)  is  evaluated  from  values  at  point  '<  at  the 

n^  iteration  level.  In  the  determination  of  v^n^  and  tan  ,  Lagrange 

cubic  interpolation  is  used.  Convergence  of  u  is  quite  rapid,  requiring 

usually  about  four  or  five  iterations  to  reach  | Au  |  <  10  Once  (u,v) 

K  K. 

are  known,  are  found  by  the  inverse  transformation  of  the  mapping 

sequence . 


Stretching  Functions 

In  the  x  -  6  plane  the  two-sided  stretching  function  of  Vinokur  is  used 
to  generate  the  grid  point  distributions  on  the  stagnation  streamline  H-E 
ind  the  airfoil-wake  centerline  E-K.-F-I.  For  segment  B-K  a  single  stretehin 
tunc  cion  is  used  whereas  for  segment  K-K-F-I,  a  sequence  of  three  stretching 


functions  is  required. 
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The  two-sided  stretching  function  for  the  normalized  variable  t  is 
given  by 


_ tanh(  ) _ 

A  sinh  A4>  +  (1  -  A  cosh  A<J> )  tanh(  ) 


(30) 


where  £  is  the  normalized  generating  variable  of  constant  step  size,  A<J>  is 
the  solution  of  the  transcendental  equation 


and 


B 


sinh  Aj> 
A<j> 


(31) 


A  =  (Sq/S^172  (32) 

B  =  (SqSj ) 1/2  (33) 


and  Sq  and  are  dimensionless  slopes  defined  as 


so-;£co)  • 


s 


1 


> 


which  control  the  clustering  at  t  =  0  and  t  =  1. 

The  reason  for  using  two-sided  stretching  functions  on  segments  E-K,  K-F 
and  F-L  is  to  provide  clustering  at  all  segment  end  points.  It  is  needed  at 
point  E  because  ot  the  rapid  drop  in  pressure  downstream  of  the  stagnation 
point,  at  point  K  to  provide  a  more  nearly  uniform  grid  distribution  on  the 
inflow  line  and  at  points  F  and  I  to  resolve  the  flow  at  the  airfoil  trailing 
edge  and  tai L  of  the  bodv  respectively.  Since  the  arc  length  step  size  at 
points  K  and  F  should  be  continuous,  not  all  of  the  parameters  are 
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independent.  If  (Sq,S^)  are  specified  on  E-K  and  is  specified  on  K-F, 
then  Sq  on  these  latter  segments  must  be  calculated  to  satisfy  continuity  of 
step  size  at  the  segment  junctions.  The  relation  is 


(so’kf  ■  (S1> 


EK 


KF 

5EK 


EK 

Vf 


(34) 


where  s  denotes  arc  length  of  the  segment  and  N  the  number  of  intervals  on 
the  segment.  A  similar  expression  holds  on  segment  F-I. 


Transfinite  Interpolation 

Normalized  pseudo-computational  variables  £  and  n  are  defined  such  that 
the  interior  of  the  quadrilateral  in  Fig.  10  in  the  u  -  v  plane  transforms 
to  the  interior  of  the  unit  square  in  the  £  -  n  plane.  The  transformation 
from  the  computational  domain  to  the  hinge-plane  domain  is  given  in  terms  of 
the  position  vector  r: 


r(5,h) 


u(C,n) 


v(?.n) 


I 


(35) 


where  0<f<l,0<q<l. 

Specifying  the  distribution  of  the  position  vector  r  and  its  normal 
derivatives  on  the  four  boundaries  in  the  C  -  h  plane  is  equivalent  to 
defining  the  grid  on  the  boundaries  in  the  hinge  plane  and  ultimately, 
the  x  -  0  plane. 
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The  transfinite  interpolation  method  used  here  is  the  extension  of 
Eriksson  [1]  as  specialized  by  Vinokur  and  Lombard  [5].  The  relation  for 
r,  using  point-of-compass  notation  for  the  boundaries,  is 

r(tn)  =  rs(c)E(n)  +  *N(5)F(n)  +  \  (c)G(n)  +  ^  (e)H(n) 


+  F(5)[?F(n)  -  r_pE( o)  -  t  F(n)  -  K  G(n)  -  i  H(n)] 


+  (n)  -  v  e(S)  -  t.  F(n)  -  r  c(n)  -  r  «(n)] 

^SW  ^NW  ^  SW  ^ nNW 


+  (n)  -  r  E(n) 

^E  oE 


r  F(n)  -  r  G(n) 

^NE  *nSE 


*  'ne 


where  E,  F,  G  and  H  are  cubic  blending  functions  given  by 


F(u)  =  u  (3  -  2u) 


G(u)  =  u( 1  -  u)' 


H(u)  =  u  (u  -  1) 


E(u)  =  l  -  F(u) 


'•*.*•*  •-  .  . 
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Equation  (36)  thus  provides  a  smooth  blending  on  the  interior  of  the  given 
distribution  of  grid  points  and  normal  derivatives  on  the  boundaries.  A 
typical  grid  in  the  hinge  plane  obtained  by  transfinite  interpolation  is 
shown  in  Fig.  11. 

The  evaluation  of  the  various  derivatives  on  the  boundaries  in  Eq.  (36) 
follows  the  prescription  given  by  Vinokur  &  Lombard  and  is  presented  in 
detail  in  Appendix  II. 

Addon  Grid 

The  foregoing  procedure  produces  a  C-grid  in  the  x  -  0  plane  in  the 

A 

region  upstream  of  the  tail  line  I-J.  Because  £  =  constant  lines  in  the 
upstream  grid  are  normal  to  I-J  (and  I-J  is  straight  as  well  as  normal  to  the 
wake  centerline),  a  downstream  grid  can  easily  be  created  which  has  continuity 
through  first  derivatives  across  I-J.  The  addon  grid  which  has  these 
characteristics  is  a  Cartesian  grid  with  the  same  0  distribution  at  I-J  as 
the  upstream  grid.  Distributing  grid  points  on  the  x-direction  downstream 
of  I-J  is  accomplished  by  a  one-sided  Vinokur  stretching  function  with  the 
parameter  Sq  determined  by  requiring  continuity  of  Ax  on  either  side  of  I-J. 

2 .4  Computational  Grid 

If  indices  i,  j  and  k  denote  the  coordinates  £,  n  and  r,  then  the 
computational  coordinates  x,  v  and  z  may  be  conveniently  defined  as 

X  =  i  -  l  ,  1  <  i  <  L 

max 

Y  =  i  -  l  ,  1  <  j  <  j  !>  .  (38) 

•'  ’  J  Jmax  i 


Z  =  k  -  1 


l  <  k  <  k 
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3 .  Results  and  Discussion 

The  afterbody-fin  grid  generation  code  is  called  TAILGRID  and  consists  of 

about  1600  FORTRAN  statements.  It  is  written  in  double  precision  arithmetric 

and  computes  in  terms  of  real  variables  only.  To  date  all  grid  generation 

has  been  done  on  a  VAX  11/782  computer  with  CPU  per  grid  point  found  to  be 
-3 

about  7  x  10  sec.  Thus  computing  a  surface  containing  1500  points  requires 
approximately  10  sec. 

The  airfoil  family  chosen  for  testing  the  grid  generation  procedure  was 
the  NACA  symmetric  four  digit  series.  The  equation  for  this  profile  is 

yp  =  -  5t(0 .2969  /7  -  0.1281x  -  0.3516x2  +  0.2843x3  -  0.1015x4)  ,  (39) 

where  x  is  the  maximum  thickness  expressed  as  a  fraction  of  the  chord.  In 

the  original  equation  for  yp ,  see  Eq.  (6.2)  of  Ref.  8,  the  coefficient  of  x 

is  given  as  0.12600  which  causes  the  airfoil  to  have  a  finite  trailing  edge 

thickness  (yTp  =  0.0021 ).  Since  the  grid  generation  procedure  requires  zero 

trailing  edge  thickness,  the  coefficient  of  x  was  modified  as  shown  in 

Eq.  (39).  Interpolation  is  used  liberally  on  the  airfoil  in  the  grid 

generation  process;  thus  an  accurate  definition  of  y„  versus  x  is  a 

F 

necessity.  Usually  100  points  on  the  airfoil  are  computed  for  this  purpose 
with  clustering  at  the  leading  edge. 

A  number  of  2-D  test  cases  in  the  x  -  0  plane  were  run  to  determine 
the  effect  of  certain  input  parameters  on  grid  quality.  These  cases  all 
consisted  of  an  NACA  0012  airfoil  at  a  cylindrical  radius  of  0.5.  The 
common  parameters  for  the  five  cases  are  given  in  Table  1  while  the 
parameters  that  vary  from  case  to  case  are  presented  in  Table  2. 
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N1 

N2 

JMAX 

DOB 

SXO 

SX2 

SX3 

SY1 

10 

20 

31 

1.0 

5.0 

2.0 

1.1 

1 .0 

Table  1.  Common  Parameters  in  2-D  Test  Cases. 


Table  2.  Variable  Parameters  in  2-D  Test  Cases. 

In  all  of  the  above  cases,  the  point  distributions  on  inflow-lower  symmetry 
line  and  outflow  boundaries  are  determined  as  described  in  the  previous  section. 
In  this  2-D  example  only  the  grid  point  distribution  downstream  of  the  trailing 
edge  is  given  by  a  geometric  progression. 

Case  1  and  2,  shown  in  Figs.  12  and  13,  are  the  same  except  that  clustering 
is  used  about  point  K  in  case  2  and  none  is  used  in  case  l.  The  orthogonality 
constraint  is  seen  to  produce  considerable  spreading  of  £  =  constant  lines 


28 


12  April  1985 
GHHtlhz 

near  the  corner  (point  C) .  Without  clustering  a  poor  boundary  point 
distribution  on  the  IVL  is  obtained  whereas  with  clustering  the  distribution 
close  to  point  C  is  improved. 

Case  3,  presented  in  Fig.  14,  is  similar  to  case  2  except  that  the 
clustering  parameter  on  the  stagnation  line  at  the  airfoil  leading  edge  is 
25  times  larger.  The  result  is  a  much  more  dense  grid  near  the  airfoil 
surface  which  would  be  useful  for  calculation  of  viscous  flow  at  high 
Reynolds  number.  In  both  cases  the  number  of  points  on  the  stage  ■  ion  line 
is  the  same  —  31. 

Cases  4  and  5  (see  Figs.  15  and  16)  show  what  happens  to  the  grid  when 
djV^  is  increased,  all  other  parameters  being  the  same  as  in  case  3. 

Point  K  is  seen  to  migrate  toward  the  airfoil  leading  edge  as  d.^^  increases 
which  results  in  a  squeezing  of  grid  lines  between  the  leading  edge  and 
point  K.  In  terms  of  grid  quality  the  optimum  value  of  d^^  appears  to  lie 
between  0.5  and  1.0  so  that  line  C-K  remains  nearly  straight.  One  way  of 
moving  the  IVL  further  upstream  without  sacrificing  grid  quality  would  be  to 
add  a  cartesian  grid  upstream  of  the  flattened  C-grid  with  the  same  9-spacing 
as  on  line  B-C  (the  IVL  for  the  C-grid). 

For  a  3-D  test  problem  the  afterbody  meridian  profile  was  represented 
by  the  following  analytic  function: 


rb(p)  =  r^F(y)  -  dbQdtan  Q^p) 


(40) 


where 
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V  =  ~ 


dfood  =  afterbody  length 


r^  =  initial  afterbody  radius  , 


and  F  and  G  are  the  cubic  blending  functions  defined  by  Eq.  (37).  The 
particular  values  chosen  for  the  afterbody  parameters  are 


r.  =  0.75  ,  d,  ,  =  2.5  ,  tan  0_  =  0.50  , 

bj  bod  T 


which  produce  a  fairly  full  profile  with  a  tail  half  angle  of  26.6  degrees 
Other  parameters  in  the  test  problem  are  given  in  Table  3. 


N1 

N2 

N3 

N4 

JMAX 

KMAX 

D1VL 

DTL 

DOB 

10 

20 

5 

10 

31 

3 

0.75 

0.5 

1.0 

RTIP 

TAU 

sxo 

SX1 

SX2 

SX3 

SY0 

SY1 

SY0 

1.0 

0.2 

5.0 

3.0 

2.0 

2.0 

2.0 

1.0 

2.0 

[tli] 


30 


12  April  1985 
GHH: lhz 


The  meridian  plane  view  of  the  test  problem  geometry,  computational  domain 
and  intermediate  surface  is  shown  in  Fig.  17.  Two  views  of  each  r  =  constant 
surface  are  presented,  the  first  from  below  and  in  front  and  the  second  from 
the  side.  These  views  are  shown  in  Figs.  18  through  23.  A  composite  side 
view  showing  the  position  of  each  surface  relative  to  the  other  is  presented 
in  Fig.  24. 

In  a  grid  stacking  procedure  each  grid  on  a  surface  is  generated  somewhat 
independently  of  the  other.  The  dependence  is  indirect  through  the  geometry 
and  not  direct  as  in  the  case  of  partial  differential  equation  grid  generation 
schemes  or  fully  3-D  algebraic  schemes.  Thus  for  3-D  grids  generated  by 
stacking  one  of  the  primary  concerns  is  with  smoothness  in  the  stacking 
direction.  In  the  present  method,  the  only  reason  that  the  grid  changes  in 
the  x  -  6  plane  from  surface  to  surface  is  that  the  airfoil  image  is 
changing.  As  r  increases,  the  airfoil  image,  according  to  Eq.  (17),  is 
shrinking  in  terms  of  maximum  thickness  approximately  as  1/r  .  Although  the 

d 

total  arc  length  of  the  airfoil  image  is  also  shrinking  slightly  as  r 
increases,  the  clustering  parameters  are  fixed  and  hence  the  airfoil  point 
distribution  on  each  surface  is  always  in  the  same  proportion.  As  the 
airfoil  image  grows  thinner,  point  K  slowly  moves  toward  the  leading  edge. 

The  trace  of  point  K  in  the  meridian  plane  is  shown  in  Fig.  17.  On  the  other 
hand,  the  distribution  of  points  on  the  stagnation  line  remains  the  same 
independent  of  r.  Of  course  close  to  the  surface  of  the  afterbody  in  the 
vicinity  of  the  tail  the  grid  shrinks  rapidly  to  reflect  the  pointed  nature 
of  the  tail  and  the  axis  singularity.  This  feature  would  exist  whether  or 
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not  stacking  were  used.  Thus  because  the  airfoil  image  is  varying  slowly 
as  r  increases  and  th..  oints  on  the  image  remain  in  the  same  proportion  of 
arc  length,  the  present  method  can  be  expected  to  produce  a  grid  of  high 
quality  in  the  stacking  direction. 

One  aspect  of  the  current  strategy  of  point  placement  on  boundaries  is  not 
entirely  satisfactory.  Although  the  circular  arc  method  of  point  placement 
in  the  hinge  plane  works  well  on  the  lower  symmetry  line,  it  leaves  something 
to  be  desired  on  the  IVL.  Coupled  with  the  singularity  at  the  corner 
(point  C)  and  the  orthogonality  requirement  at  boundaries,  clustering  at 
point  K  was  found  to  be  necessary  to  achieve  a  reasonable  point  spacing 
near  point  C  on  the  IVL.  This  clustering  would  probably  not  be  necessary 
if  a  different  strategy  were  used  to  locate  the  points  on  the  IVL.  One 
possibility  would  be  to  space  them  in  the  x  -  9  plane  in  the  same  proportion 
of  arc  length  as  along  the  airfoil  between  the  leading  edge  and  point  K. 
Downstream  of  point  K  the  strategy  would  remain  as  before. 
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Figure  22.  Front  View  of  Outer  Surface 
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Figure  24.  Composite  View  of  Afterbody,  Intermediate  and  Outer  Surface  Grids. 
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(2)  Transformation  from  £  -  n  to  x  -  0  plane. 


where 


C  =  y  C  +  [j  (p  +  p)]1/2 
H=|ii+  [|(p  -p)] 


y 


(p2  +  4q2) 


1/2 

9 


<i  =  |  In 

and 

x  =  Jtn(cosh  £  ~  cos  n) 


(AI.9) 

(AI.10) 

(AI.ll) 

(A1.12) 

(AI.13) 

(AI.14) 
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cos 


-1 ( 1  -  cosh  £  cos  n 
cosh  £  -  cos  n 


( 


) 


(A1.15) 


In  the  preceding  relations  £q 


is  computed  from 


£ 


0 


sinh 


[ 


4  exp(-a) 

1  -  exp(-a) 


] 


(AI.16) 


where 


a  =  4  f  d_  +  dT1„  1 


(41.17) 


59 


12  April  1985 
GHH: lhz 


In  the  determination  of  dg,  the  same  procedure  is  followed  here  as  in 
Ref.  6  where  near  the  airfoil  leading  edge  in  the  x  -  9  plane  the  most  nearly 
orthogonal  grid  is  sought.  Such  a  system  is  obtained  when  the  leading  edge 
maps  into  an  n  =  constant  line.  The  resulting  expression  for  dg  is  found  to 
be 

1  1  +  4pn 

ds  ■  4 

where  is  the  radius  of  curvature  of  the  airfoil  leading  edge  in  the  x  -  0 
plane.  Equation  (AI.18)  differs  from  its  counterpart  in  Ref.  6,  Eq .  (100), 
because  a  uniform  stretching  of  x  and  6  is  used  here. 
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Appendix  II.  Evaluation  of  Boundary  Derivatives  for  Transfinite  Interpolation 


Normalized  coordinates  |  and  n  are  determined  in  the  hinge  plane  in 
terms  of  the  normalized  arc  length  along  the  south  and  east  boundaries 
(lines  E-K-F-I  and  B-E  in  Fig.  10).  If  s^  and  S£  denote  running  arc  length 
along  the  south  and  east  boundaries,  then 

K). 

t  =  - -  ,  ( All .  1 ) 

1  SE-I 

’  (AII-2) 
where  the  running  arc  lengths  are  determined  using  the  chord  approximation. 

->  -f  A  A 

If  t^  and  t^  are  unit  tangent  vectors  along  the  £  and  n  coordinate  curves, 

then 

3s. 

r_  =  -  t.  ,  (All. 3) 

4  3C 

and 


(All. 4) 


Under  the  assumption  that  extended  orthogonality  holds  at  the  corners  the 
cross-derivative  is  given  by 


r 

A  A 


3s , 


A 

35 


3h 


* 


2 


(All. 5) 
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A  A 

where  and  are  the  curvatures  of  the  E,  and  h  coordinate  curves 
respectively. 

The  south-north  and  west-east  boundaries  are  assumed  always  to  be  oriented 
as  shown  in  Fig.  10.  Thus  on  the  south  boundary  we  take  v^  =  Vg(u)  from 
which  the  unit  tangent  vector  follows  as 


t.  =  t,  e  +  t  e 
1  lu  u  lv  v 


(All. 6) 


where 


1 


lu  .  1/2 

(i  *  <2) 


(All. 7) 


'lv  ,  1/2 

(l  +  v’2) 


(All. 8) 


and  e^  and  e^  are  unit  vectors  in  the  u  and  v  directions  respectively.  Then 
the  unit  normal  vector  to  the  south  boundary  may  be  written  as 


t_  =  -  t.  e  +  t„  e 
2  lv  u  2u  v 


(All. 9) 


And  finally,  the  curvature  on  the  south  boundary  is 


Cl  = 


17T 


(All. 10) 


is 


On  the  east  boundary  we  take  u  =  u  (v)  and  hence  the  unit  tangent  vector 

E  E 


where 
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t0  =  t„  e  +  t0  e 
2  2u  u  2v  v 


(All. 11) 


(All. 12) 


Then  the  unit  normal  vector  is 

t  =  t.  e  -  t„  e  , 
1  2v  u  2u  v 

and  the  curvature  is 


f  r 


(All. 13) 


(All. 14) 


(All. 15) 


On  the  south  boundary,  from  Eq.  (AII.l) 

9s 

( — -)  =  s  .  (All. 16) 

3?  S 

The  derivative  Ds0/9n  on  the  south  boundary  is  known  only  at  the  end  points 
and  hence  must  he  determined  between  points  E  and  I  by  interpolation. 


FoLLowing  Ref.  5,  we  use  a  cubic  blending  function  approximation,  viz. 
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,9s2 
l  „  > 
3n 
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3s„  3s  3  s 

(  ~)  e(5)  +  (-~)  f(5)  +  (-^-4)  G(c) 
3n  sv}  3n  se  H3t)  sw 


3S 

♦  (~rr)  »(« 

3?3p  spr 


(All. 17) 


From  the  extended  orthogonality  relations,  the  cross-derivative  is  found  to 
be 


A  '  A 

3?  3n 


3s^  3s2 


3?  3n 


On  the  east  boundary,  from  Eq.  (All. 2) 


(All. 18) 


hr)  -Ve  •  (AII'19) 

8n  E 


and  analogous  to  Eq.  (All. 17),  the  cubic  blending  function  formula  for 
3s^/3?  is 

2  2 

3s,  3s,  3s,  3s,  3s, 

(-r)  -  (-r1)  hS)  +  hr)  fCS)  *  h4)  c(S)  .  (_i)  ,.(S)  .  (aii.20) 

95  e  3?  SE  3?  NE  3h3C  SE  3ri3?  NE 


where  from  extended  orthogonality 


II 

7) 

fX) 

,i!i 

A  ' 

3? 

3? 

A 

3n 

(All. 21) 


To  maintain  a  reasonable  spacing  of  (  =  constant  grid  lines  in  the 
physical  Diane,  especially  near  the  airfoil  surface,  the  normal  spacing 
in  the  x  -  0  plane  on  the  stagnation  and  outflow  lines  (east  and  west 
boundaries  in  the  hinge  plane)  is  taken  to  be  the  same.  As  a  result  (p 
does  not  obey  a  relation  similar  to  Eq .  (AIT. 2)  and  hence  3sn/3p  must  be 
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determined  numerically  point-by-point.  As  for  Ss^/SC,  an  interpolation 
relation  similar  to  Eq.  (All. 20)  is  used  with  appropriate  changes  in 
notation. 

On  the  north  boundary  the  points  are  positioned  in  the  hinge  plane  by 
using  a  circular  arc  approximation  similar  to  that  used  in  locating  point  K. 
The  iteration  formula  for  u^,  similar  to  Eqs.  (28)  and  (29)  is 

K)rn  ■  + •  <aii-22> 

where 

US  ~  UN  +  tVS  ~  VN)tan[j  C*s  + 

1  +  tan  <^N  tan[y  (<f>s  +  <j>N)] 

.  * 

The  derivative  3s^/3C  must  be  determined  numerically  at  each  point  and 
3s2/3h  is  computed  from  an  interpolation  formula  analogous  to  Eq.  (All. 17). 

In  the  determination  of  the  various  interpolated  values  of  u  and  v  and 
their  derivatives  required  in  the  foregoing  equations,  Lagrange  cubic 
interpolation  is  used.  The  derivative  expressions  are  determined  by 
differentiation  of  the  Lagrange  polynomials. 


►J] 


(n) 


(All. 23) 
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